home *** CD-ROM | disk | FTP | other *** search
/ Shareware Grab Bag / Shareware Grab Bag.iso / 090 / pctj8505.arc / MARQFIT.BAS < prev    next >
BASIC Source File  |  1986-09-14  |  34KB  |  669 lines
  1. 1 ' PROGRAM MARQFIT.BAS -- a non-linear least squares fitting
  2. 2 ' routine using the marquardt algorithm
  3. 5 ' (c) 1985 W. Schreiner, M. Kramer, S. Krischer, & Y. Langsam
  4. 20  '-------------- initialize program settings ---------------
  5. 30  REV = 2.1 : PI = 3.14159265#
  6. 40  ON ERROR GOTO 19000 : RETCOD = 0 :KEY OFF
  7. 60  DEF SEG=64:POKE 23,(PEEK(23) OR 64): DEF SEG 'set caps lock
  8. 65  ' The values MXVAR%, MXOBS% and MDI% can be changed-
  9. 70  ' be sure to change the DIM sizes correspondingly
  10. 80  MXVAR% = 30:MXOBS% = 500 : MDI% = MXVAR% * (MXVAR% + 1) / 2
  11. 90  DIM X(500),DTA(500),WGT(500),FXP(500),FYP(500),XP(500)
  12. 100 DIM YP(500),C(465),A(465),B(30),DPARM(30),PARM(30),                                 KFIX(30),DERIV(30),G(30),GRAD(30)
  13. 110 DEF FNLGT(X) = LOG(X)/LOG(10)
  14. 120 ON KEY(1) GOSUB 10200 : ON KEY(2) GOSUB 10300
  15. 130 ON KEY(3) GOSUB 10310 : ON KEY(4) GOSUB 10400
  16. 140 KEY(1) ON : KEY(2) ON : KEY(3) ON : KEY(4) ON
  17. 160   FG=2:BG=0:FG1=14:BG1=0:FG2=4:BG2=0 ' color attributes
  18. 165 FOR I = 1 TO MXOBS% : WGT(I) = 1 : NEXT
  19. 190 ' ---------------plot routine initialization
  20. 200 PLOTSET% = 0 ' plot definition flag
  21. 210 XTOTAL.PIX = 639: YTOTAL.PIX = 199
  22. 220 ' (for lo_res use XTOTAL.PIX = 319)
  23. 230 ' XMIN.PIX is lower left hand xmin in pixel coords;
  24. 240 ' YMIN.PIX is lower left hand ymin in pixel coords;
  25. 250 ' XMIN and YMIN are the same but in USER coordinates
  26. 260 XMIN.PIX=CINT(.11*XTOTAL.PIX): XMAX.PIX=.94*XTOTAL.PIX
  27. 270 YMAX.PIX=CINT(.07*YTOTAL.PIX)-2: YMIN.PIX=YTOTAL.PIX-12
  28. 280 DELTAX.PIX = XMAX.PIX-XMIN.PIX+1
  29. 290 DELTAY.PIX = YMIN.PIX-YMAX.PIX+1
  30. 310 ' functions to convert X & Y in user coords --> pixel coords
  31. 320 DEF FNSX(N)=XMIN.PIX + CINT((N-XMIN)*DELTAX.PIX/(XMAX-XMIN))
  32. 330 DEF FNSY(N)=YMIN.PIX - CINT((N-YMIN)*DELTAY.PIX/(YMAX-YMIN))
  33. 335 ' ------------------------start up--------------------------
  34. 340 IPFLG = 0    ' pause-in-fitting flag
  35. 350 ISVFL =  - 1 ' data saved flag
  36. 360 CLS : LOCATE 10 : COLOR FG2,BG2 :                                               PRINT TAB(15) "MARQUARDT LEAST SQUARES - REV ";REV:
  37. 365 COLOR FG,BG : FOR I = 1 TO 1700: NEXT:LOCATE 20
  38. 370 INPUT "USE DATA ON DRIVE (ENTER LETTER):";DRV$:                                 DRV$ =LEFT$(DRV$,1) : IF DRV$ <> "" THEN DRV$=DRV$+":"
  39. 380  GOTO 420
  40. 390 ' -----------------------   menu  --------------------------
  41. 400 PRINT "PRESS ANY KEY TO CONTINUE":WHILE INKEY$ = "":WEND
  42. 410 RETCOD = 0 ' reset error flag
  43. 420 F1%=0:SCREEN 0 :CLS: GOSUB 10000: COLOR FG2,BG: LOCATE 1,20
  44. 425 PRINT "MENU OPTIONS:": COLOR FG,BG : PRINT:COLOR FG1,BG1
  45. 430 PRINT TAB(6) "GENERAL:";:COLOR FG,BG
  46. 435 PRINT TAB(21) "1- ENTER TITLE " 'line  800
  47. 440 PRINT  TAB(21) "2- PRINTER ";   'line  850
  48. 450 IF PFLG% = 0 THEN PRINT "ON/";:COLOR FG1,BG1:PRINT  "OFF":                      COLOR FG,BG ELSE COLOR FG1,BG1: PRINT "ON";: COLOR FG,BG :                      PRINT  "/OFF"
  49. 460  PRINT  TAB(21) "3- SOLVE "      'line  900
  50. 470  PRINT  TAB(21) "4- PLOT "       'line  15000
  51. 480  PRINT  TAB(21) "5- QUIT "       'line  1300
  52. 490  PRINT  : COLOR FG1,BG1:PRINT  "     ENTER DATA:";:                              COLOR FG,BG:PRINT TAB(21) "6- MANUAL" 'line 1400
  53. 500  PRINT  TAB(21) "7- UREAD"       'line 1800
  54. 510  PRINT
  55. 520  COLOR FG1,BG1:PRINT  "     MODIFY DATA:";:COLOR FG,BG:                          PRINT TAB(21) "8- EDIT"         'line 2000
  56. 530  PRINT  TAB(21) "9- SCALE"       'line  2200
  57. 540  PRINT  TAB(20) "10- ZERO"       'line  2400
  58. 550  PRINT:COLOR FG1,BG1:PRINT "     PARAMETERS:";:COLOR FG,BG:                      PRINT TAB(20) "11- ENTER/REVIEW/CHANGE"     'line 2500
  59. 570  PRINT  "                   12- FIX"         'line  2800
  60. 580  PRINT  "                   13- FREE"        'line  3000
  61. 590  PRINT
  62. 600  COLOR FG1,BG1:PRINT  "     DATA & PARAM: ";:COLOR FG,BG:                        PRINT "14- LIST "     ' line 3200
  63. 610  PRINT  "                   15- SAVE ON DISK" 'line 3400
  64. 620  PRINT  "                   16- READ FROM DISK" 'line 3600
  65. 630  COLOR 15,BG: INPUT; "ENTER CHOICE";KMND:COLOR FG,BG:PRINT
  66. 650  ON KMND+1 GOTO 420,800,850,900,15000,1300,1400,1800,2000,                       2200,2400,2500,2800,3000,3200,3400,3600
  67. 660  PRINT  "*** ERROR *** ";KMND;" INVALID COMMAND": GOTO 400
  68. 800 'enter a title for documentation ---------------------------
  69. 810 INPUT "ENTER TITLE: ";TITL$:LOCATE 25,1:PRINT TITL$;SPC(11)
  70. 820  GOTO 420
  71. 850 'toggle printer on/off -------------------------------------
  72. 860 PFLG% = 1 - PFLG%: GOTO 420
  73. 900 'solve the mrqdt non-linear lst sq fit problem -------------
  74. 920 ITER% = 0 : IF IPFLG <  > 0 THEN ITER% = NITER%
  75. 930 IPFLG = 0:ISVFL = 0 : CLS : GOSUB 10000
  76. 940 INPUT "HOW MANY ITERATIONS?  ";MXITER%
  77. 950 IF MXITER% < 0 THEN MXITER% = 0
  78. 960 GOSUB 6000 ' go fit !
  79. 970 IPFLG =  - 1 : NITER% = NITER% + ITER%
  80. 990 PRINT:PRINT: PRINT TITL$ :PRINT NITER%;" ITERATIONS": PRINT
  81. 1000 IF PFLG% = 1 THEN LPRINT:LPRINT TITL$                                            :LPRINT NITER%;" iterations "
  82. 1010 GOSUB 9000
  83. 1020 IF FLG%=1 THEN PRINT:PRINT "CHOLESKY NEGATIVE DIAGONAL--";                            "UNABLE TO  SOLVE WITH SUPPLIED INITIAL PARAMETERS"
  84. 1030 PRINT  : PRINT  "PRESS ANY KEY FOR FITTING STATISTICS"
  85. 1040 WHILE INKEY$ = "" : WEND
  86. 1050 GOSUB 1100 : PRINT : GOTO 400 'print statistics and return
  87. 1070 'print fit statistics -------------------------------------
  88. 1100 PRINT  : PRINT  "SOME FITTING STATISTICS:"
  89. 1110 PRINT "SIGMA= ";SIG;" R= ";R;: PRINT TAB( 41);                                                                    "WGT'D SUM OF SQ. = ";WSS
  90. 1120 CHI2 =  INT (10 * CHI2) / 10
  91. 1130 PRINT  "CHI SQUARED= ";CHI2;" / ";NF;" DEG OF FREEDOM"
  92. 1140 PRINT  "# OF CALLS TO SUM OF SQ=";NSSC;:                                             PRINT TAB( 41);"# OF DERIV CALLS=";NDC:                                         PRINT  "# INC IN LAMBDA=";INCR;" LAMBDA= ";LAMBDA
  93. 1150 IF PFLG% <> 1 THEN RETURN
  94. 1160 LPRINT : LPRINT  "SOME FITTING STATISTICS:"
  95. 1170 LPRINT  "SIGMA= ";SIG;" R= ";R;:                                                 LPRINT TAB( 41);"WGT'D SUM OF SQ. = ";WSS
  96. 1180 LPRINT  "CHI SQUARED= ";CHI2;" / ";NF;" DEG OF FREEDOM"
  97. 1190 LPRINT "# OF CALLS TO SUM OF SQ=";NSSC;:                                        LPRINT   TAB( 41);"# OF DERIV CALLS=";NDC:                                      LPRINT  "# INC IN LAMBDA=";INCR;" LAMBDA= ";LAMBDA
  98. 1200 RETURN
  99. 1300 'quit -----------------------------------------------------
  100. 1310  IF ISVFL <  > 0 THEN 1360
  101. 1320  INPUT "PRESENT DATA NOT SAVED. SAVE IT?? (Y/N)";T1$
  102. 1330  IF  LEFT$ (T1$,1) = "Y" THEN 3400 'go save
  103. 1340  IF  LEFT$ (T1$,1) <> "N" THEN 420
  104. 1360  DEF SEG=64:POKE 23,(PEEK(23) AND 191):DEF SEG 'clear caps
  105. 1370  CLS: END
  106. 1400  'manual data entry ---------------------------------------
  107. 1410  PRINT  "MANUAL DATA ENTRY - ": PRINT  "   (, TO EXIT)"
  108. 1412 ' set flags to indicate data points that currently exist
  109. 1415 FOR K=1 TO NOBS%:XP(K)=1:NEXT :                                                  FOR K = NOBS%+1 TO MXOBS%:XP(K)=0:NEXT
  110. 1420 PRINT  "POINT    X,MEASUREMENT"
  111. 1425 PRINT  "===== ====,==========="
  112. 1430 FOR I = 1 TO MXOBS%
  113. 1440    PRINT  I;: INPUT " ";T1$,T2$
  114. 1450    IF T1$ = "" THEN NOBS% = I - 1: GOTO 1480
  115. 1460    X(I) =  VAL (T1$):DTA(I) =  VAL (T2$) :XP(I) = 1
  116. 1470 NEXT
  117. 1480 INPUT "WEIGHTS-ENTER 'NO','STAT' OR 'EXPLICIT'?";T1$
  118. 1490 IF T1$ = "NO" THEN FOR I=1 TO NOBS%:WGT(I)=1:NEXT:GOTO 1580
  119. 1500 IF T1$ = "STAT" THEN                                                              FOR I = 1 TO NOBS%: WGT(I)=SQR (DTA(I)): NEXT : GOTO 1580
  120. 1510 IF T1$ <  > "EXPLICIT" GOTO 1480
  121. 1520 PRINT  : PRINT  "POINT     X      Y      WEIGHT":                                        PRINT  "=====   =====  =====    ======"
  122. 1530 FOR I = 1 TO NOBS%
  123. 1540     PRINT  I;"     ";X(I);"    ";DTA(I);"     ";WGT(I);
  124. 1550     INPUT "New Weight=";WGT$:IF WGT$="" THEN 1570
  125. 1560     WGT(I) = VAL(WGT$)
  126. 1570 NEXT
  127. 1580 NOBS% = 0 : PRINT  "NOW ORDERING & REVIEWING DATA"
  128. 1600 FOR I = 1 TO MXOBS%
  129. 1610    IF XP(I) = 0 GOTO 1650  'zero it if pt is to be deleted
  130. 1620    NOBS% = NOBS% + 1
  131. 1630    X(NOBS%) = X(I):DTA(NOBS%) = DTA(I):WGT(NOBS%) = WGT(I)
  132. 1640    IF I <  > NOBS% THEN XP(I) = 0
  133. 1650 NEXT
  134. 1680 IPFLG = 0:ISVFL = 0 : GOTO 420
  135. 1800 'user defined read routine --------------------------------
  136. 1810 PRINT "NO USER DEFINED ROUTINE IMPLEMENTED"
  137. 1820 GOTO 400
  138. 2000 'edit data ------------------------------------------------
  139. 2005 FOR K=1 TO NOBS%:XP(K)=1:NEXT
  140. 2010 FOR K = NOBS%+1 TO MXOBS%:XP(K)=0:NEXT
  141. 2015 PRINT "ENTER 'W' TO CHANGE ONLY WEIGHTS,";:                                                      INPUT "'D' TO CHANGE DATA & WEIGHTS";T1$
  142. 2020 IF T1$ = "W" GOTO 1480
  143. 2030 PRINT "ENTER DATA PT RANGE FOR EDITING (N1,N2)"
  144. 2040 INPUT "(, TO EXIT)";T1$,T2$
  145. 2050 IF T1$ = "" GOTO 1580
  146. 2060 I =  VAL (T1$):J =  VAL (T2$)
  147. 2070 IF I<1 OR J>MXOBS% THEN PRINT "INVALID RANGE":GOTO 2030
  148. 2080 FOR K = I TO J
  149. 2090     PRINT  "CURRENT VALUE OF X,DTA,WGT FOR PT #";K;" =":                               PRINT X(K);"  ";DTA(K);"  ";WGT(K)
  150. 2100     PRINT  "ENTER NEW VALUE FOR X(";K;")"
  151. 2105     PRINT "(ENTER 'D' TO DELETE PT, OR 'enter' ";:                                  INPUT "TO LEAVE UNCHANGED)";T1$
  152. 2110     IF T1$ = "D" THEN XP(K) = 0: GOTO 2150
  153. 2120     IF T1$ = "" THEN 2150
  154. 2130     X(K) =  VAL (T1$) : XP(K) = 1
  155. 2140     INPUT "ENTER NEW VALUES FOR DTA,WGT ";DTA(K),WGT(K)
  156. 2150 NEXT
  157. 2160 GOTO 2030
  158. 2200 'scale x,dta,wgt --------------------------------------
  159. 2210 INPUT "ENTER X-COORDINATE SCALE FACTOR";XSCALE
  160. 2220 INPUT "ENTER DATA POINT SCALE FACTOR";DSCALE
  161. 2230 INPUT "ENTER WEIGHT SCALE FACTOR";WSCALE
  162. 2240 FOR I = 1 TO NOBS%
  163. 2250 X(I)=X(I)*XSCALE:DTA(I)=DTA(I)*DSCALE:WGT(I)=WGT(I)*WSCALE
  164. 2280 NEXT
  165. 2290 IPFLG = 0:ISVFL = 0 : GOTO 420
  166. 2400 'zero data & parameters -----------------------------------
  167. 2410 FOR I=1 TO MXOBS%:X(I)=0:DTA(I)=0:WGT(I)=0:NEXT:NOBS% = 0
  168. 2420 FOR I=1 TO MXVAR%:PARM(I)=0:DPARM(I) = 0: NEXT :NPRM% = 0
  169. 2430 IPFLG = 0:ISVFL =  - 1 : GOTO 420
  170. 2500 'enter/review/change param --------------------------------
  171. 2510 CLS:GOSUB 9000 : IF NPRM% < 1 THEN 2580
  172. 2530 PRINT  "CHANGE PARAMETERS? ENTER:   'A' TO CHANGE ALL":                        PRINT TAB(29)"'S' TO SELECTIVELY CHANGE (or add parameters)"
  173. 2540 PRINT TAB(29) "press RETURN to leave unchanged";:INPUT T1$
  174. 2550 IF T1$ = "" GOTO 420
  175. 2560 IF T1$ = "S" THEN 2700
  176. 2570 IF T1$ <  > "A" THEN 2530
  177. 2580 PRINT "ENTER PARAMETERS ONE AT A TIME.":                                            PRINT "   (NULL TO EXIT)"
  178. 2590 FOR I = 1 TO MXVAR%
  179. 2600     PRINT  I;"- "; : INPUT T$ : PRINT
  180. 2610     IF T$ = "" THEN NPRM% = I - 1: GOTO 2660
  181. 2620     IF  LEFT$ (T$,1) = "." OR LEFT$ (T$,1) = "-" THEN 2640
  182. 2630     IF  LEFT$ (T$,1) < "0" OR  LEFT$ (T$,1) > "9" THEN                                 PRINT T$ "INVALID, RETYPE": GOTO 2600
  183. 2640     PARM(I) =  VAL (T$)
  184. 2650 NEXT
  185. 2660 IF NPRM% <  = 0 GOTO 420
  186. 2670 FOR I = 1 TO NPRM%:DPARM(I) = 0: NEXT
  187. 2680 IPFLG = 0 : ISVFL = 0 : GOTO 420
  188. 2700 '----- change some of the param
  189. 2710 PRINT "ENTER PARM#, AND VALUE": PRINT SPC(10)"(, TO EXIT)"
  190. 2720 INPUT T1$,T2$: IF T1$ = "" THEN  GOSUB 9000: GOTO 400
  191. 2730 IF VAL(T1$) = NPRM%+1 THEN NPRM%=NPRM%+1
  192. 2740 IF  VAL (T1$) < 1 OR  VAL (T1$) > NPRM% THEN                                        PRINT "INVALID PARM#. RANGE IS 1-";NPRM%+1: GOTO 2720
  193. 2750 IF  LEFT$ (T2$,1) = "." OR LEFT$ (T2$,1) = "-" THEN 2770
  194. 2760 IF  LEFT$ (T2$,1) < "0" OR  LEFT$ (T2$,1) > "9" THEN                                PRINT T2$ "INVALID, RETYPE": GOTO 2720
  195. 2770 PARM( VAL (T1$)) =  VAL (T2$)
  196. 2780  GOTO 2720
  197. 2800 'fix some parameters -------------------------------------
  198. 2810 GOSUB 9000
  199. 2820 IF NPRM% < 1 THEN 400
  200. 2830 INPUT "ENTER PARM# TO BE FIXED (NULL TO EXIT)";I
  201. 2840 IF I = 0 GOTO 2870
  202. 2850 IF I<1 OR I>NPRM% THEN PRINT "INVALID PARM#. RANGE IS 1-"                          ;NPRM%: GOTO 2800
  203. 2860 KFIX(I) = 1 : F1% = 1 : GOTO 2830
  204. 2870 IF F1% = 0 GOTO 400 'NOTHING CHANGED
  205. 2880 GOSUB 9000 : IPFLG = 0:ISVFL = 0 : GOTO 400
  206. 3000 'free some previously fixed param.-------------------------
  207. 3010 GOSUB 9000 : IF NPRM% < 1 THEN 400
  208. 3030 INPUT "ENTER PARM# TO BE FREED (NULL TO EXIT)";I
  209. 3040 IF I = 0 GOTO 3070
  210. 3050 IF I < 1 OR I > NPRM% THEN                                                          PRINT  "INVALID PARM#. RANGE IS 1-";NPRM%: GOTO 3000
  211. 3060 KFIX(I) = 0 : F1% = 1 : GOTO 3030
  212. 3070 IF F1% = 0 GOTO 400
  213. 3080 GOSUB 9000 : IPFLG = 0:ISVFL = 0 : GOTO 400
  214. 3200 'print data and param values ------------------------------
  215. 3210 CLS:GOSUB 10000:PRINT "MARQUARDT LEAST SQUARES - REV ";REV
  216. 3220 PRINT  TITL$: PRINT  NITER%;" ITERATIONS": PRINT
  217. 3230 IF PFLG%=1 THEN LPRINT TITL$:LPRINT NITER%;" ITERATIONS"
  218. 3240 GOSUB 9000 : GOSUB 1100 'print parameters and statistics
  219. 3260 PRINT:PRINT NOBS%;" DATA POINTS BEING FITTED":                                    PRINT SPC(18)"(X,DTA,WGT,FMM)"
  220. 3270 IF PFLG%=1 THEN LPRINT:LPRINT NOBS%;                                           " DATA POINTS BEING FITTED":LPRINT SPC(18) "(X,DTA,WGT,FMM)"
  221. 3280 IF NOBS% <  = 0 GOTO 3360
  222. 3290 PFORM$ = "##.##^^^^   "
  223. 3300 FOR I = 1 TO NOBS%
  224. 3310 IF BRKFLG% THEN BRKFLG% = 0:PRINT "BREAK IN LISTING":                              GOTO 400
  225. 3320 GOSUB 20000:FMM = FUNCTN - DTA(I)
  226. 3330 PRINT USING PFORM$;X(I),DTA(I),WGT(I),FMM
  227. 3340 IF PFLG%=1 THEN LPRINT USING PFORM$;X(I),DTA(I),WGT(I),FMM
  228. 3350 NEXT
  229. 3360 GOTO 400
  230. 3400 'save data  & parameters ---------------------------------
  231. 3410 PRINT  "THESE ARE YOUR CURRENT FILES" : FILES DRV$ + "*.*"
  232. 3430 PRINT  : INPUT "SAVE DATA AS DISK FILE NAMED ";T1$
  233. 3435 IF MID$(T1$,2,1) <> ":" THEN T1$ = DRV$ + T1$
  234. 3440 OPEN "O",2,T1$ : IF RETCOD <> 0 THEN CLOSE 2 : GOTO 400
  235. 3460 WRITE #2, TITL$ : WRITE #2, NOBS%
  236. 3470 FOR I = 1 TO NOBS%: WRITE #2, X(I),DTA(I),WGT(I): NEXT
  237. 3480 WRITE #2, NPRM%
  238. 3485 FOR I = 1 TO NPRM%:WRITE #2, PARM(I),DPARM(I),KFIX(I):NEXT
  239. 3490 CLOSE 2: IPFLG = 0:ISVFL =  - 1 : GOTO 420
  240. 3600 'read data & param from disk -----------------------------
  241. 3610 PRINT "THESE ARE YOUR CURRENT FILES" : FILES DRV$ + "*.*"
  242. 3630 PRINT:INPUT "READ DATA & PARAM FROM DISK FILE NAMED ";T1$
  243. 3635 IF MID$(T1$,2,1) <> ":" THEN T1$ = DRV$ + T1$
  244. 3640 OPEN "I",3,T1$ : IF RETCOD <> 0 THEN CLOSE 3 : GOTO 400
  245. 3660 INPUT #3, TITL$ : INPUT #3, NOBS%
  246. 3670 FOR I = 1 TO NOBS%: INPUT #3, X(I),DTA(I),WGT(I): NEXT
  247. 3680 INPUT #3, NPRM%
  248. 3685 FOR I = 1 TO NPRM%:INPUT #3, PARM(I),DPARM(I),KFIX(I):NEXT
  249. 3690 CLOSE 3
  250. 3700 IPFLG = 0:ISVFL =  - 1:GOSUB 9000:PRINT
  251. 3705 LOCATE 25,1:PRINT TITL$;:LOCATE 21 : PRINT  : GOTO 400
  252. 4000 'subroutine to compute chi-squared ------------------------
  253. 4010 CHI2 = 0
  254. 4020 FOR I = 1 TO NOBS%
  255. 4030    IF WGT(I) = 0 THEN 4060
  256. 4040    GOSUB 20000
  257. 4050    CHI2 = CHI2 + ((FUNCTN - DTA(I)) / WGT(I)) ^ 2
  258. 4060 NEXT
  259. 4070 RETURN
  260. 4500 'subroutine obtains CHOLESKY backward solution of matrix A
  261. 4510 G(1) = G(1) / A(1)
  262. 4520 IF NFIT% <  = 1 THEN 4630
  263. 4530 L = 1
  264. 4540 FOR I = 2 TO NFIT%
  265. 4550     K = I - 1
  266. 4560     FOR J = 1 TO K
  267. 4570        L = L + 1 : G(I) = G(I) - A(L) * G(J)
  268. 4590     NEXT
  269. 4600     L = L + 1 : G(I) = G(I) / A(L)
  270. 4620 NEXT
  271. 4630 MDI% = NFIT% * (NFIT% + 1) / 2
  272. 4640 G(NFIT%) = G(NFIT%) / A(MDI%)
  273. 4650 IF NFIT% <  = 1 THEN  RETURN
  274. 4660 FOR K1 = 2 TO NFIT%
  275. 4670     I = NFIT% + 2 - K1
  276. 4680     K = I - 1 : L = I * K / 2
  277. 4700     FOR J = 1 TO K
  278. 4710        G(J) = G(J) - G(I) * A(L + J)
  279. 4720     NEXT
  280. 4730     G(K) = G(K) / A(L)
  281. 4740 NEXT
  282. 4750 RETURN
  283. 5000 'subroutine to perform CHOLESKI decomposition of matrix a
  284. 5010 FLG% = 0
  285. 5020 FOR J = 1 TO NFIT%
  286. 5030     L = J * (J + 1) / 2
  287. 5040     IF J <  = 1 THEN 5120
  288. 5050     FOR I = J TO NFIT%
  289. 5060        K1 = I * (I - 1) / 2 + J  : F = A(K1) : K2 = J - 1
  290. 5090        FOR K = 1 TO K2:F = F - A(K1 - K) * A(L - K): NEXT
  291. 5100        A(K1) = F
  292. 5110     NEXT
  293. 5120     IF A(L) > 0 THEN 5160
  294. 5130     FLG% = 1 : PRINT "CHOLSKI-NEG DIAG J,L,A(L)= ";J,L,A(L)
  295. 5150     RETURN
  296. 5160     F =  SQR (A(L))
  297. 5170     FOR I = J TO NFIT%
  298. 5180        K2 = I * (I - 1) / 2 + J :  A(K2) = A(K2) / F
  299. 5200     NEXT
  300. 5210 NEXT
  301. 5220 RETURN
  302. 5500 'subroutine to calculate the weighted sum of squares ------
  303. 5510 WSS = 0 : NSSC = NSSC + 1
  304. 5520 FOR I = 1 TO NOBS%
  305. 5530     GOSUB 20000
  306. 5540     WSS = WSS + (WGT(I) * (FUNCTN - DTA(I))) ^ 2
  307. 5550 NEXT
  308. 5570 RETURN
  309. 5700 'subroutine to calculate the sum of squares ---------------
  310. 5710 SS = 0 : NSSC = NSSC + 1
  311. 5720 FOR I = 1 TO NOBS%
  312. 5730     GOSUB 20000 : SS = SS + (FUNCTN - DTA(I)) ^ 2
  313. 5750 NEXT
  314. 5770 RETURN
  315. 5900 'subroutine to calculate the sum of DTA(I)^2 --------------
  316. 5910 R = 0
  317. 5920 FOR I = 1 TO NOBS% : R = R + DTA(I) ^ 2 : NEXT
  318. 5950 RETURN
  319. 5990 'subroutine to do the main mathematics of the fitting -----
  320. 6000 NFIT% = 0:NX% = 0
  321. 6010 IF NPRM% < 1 THEN RETURN 2500
  322. 6020 FOR I = 1 TO NPRM%
  323. 6030     IF KFIX(I) = 0 GOTO 6070
  324. 6040     NX% = NX% + 1
  325. 6050     DPARM(NX%) = 0!
  326. 6060     GOTO 6090
  327. 6070     NFIT% = NFIT% + 1
  328. 6080     B(NFIT%) = PARM(I)
  329. 6090 NEXT
  330. 6100 IF NFIT% <= 0 THEN NITER% = ITER% : RETURN 420
  331. 6140 'perform MARQUARDT fitting --------------------------------
  332. 6150 IF NOBS% >= NFIT% THEN 6180
  333. 6160 PRINT  "# OF DATA PTS (NOBS%= ";NOBS%;")";" IS LESS THAN";                             TAB( 41);"# OF PARAM. TO BE FIT (NFIT%= ";NFIT%;")"
  334. 6165 PRINT  "MODEL IS UNDETERMINED"; CHR$ (7)
  335. 6170 NITER% = ITER% : RETURN 400
  336. 6180 PRCSN = 2.5E-07 : LAMBDA = .01 : LMIN = 1E+20 : PHI = 1
  337. 6190 MDI% = NFIT% * (NFIT% + 1) / 2
  338. 6200 NSSC = 0 : NDC = 0 : NITER% = 0 : INCR = 0
  339. 6210 GOSUB 5500 ' Choleski decomposition
  340. 6220 PRINT "WGT'D SSQ = ";WSS
  341. 6225 IF PFLG%=1 THEN LPRINT "WGT'D SSQ = ";WSS
  342. 6230 GOTO 6280
  343. 6240 'NEXT ITERATION -------------------------------------------
  344. 6250 NITER% = NITER% + 1
  345. 6260 IF (NITER% > 4 AND LAMBDA < LMIN) THEN LMIN = LAMBDA
  346. 6270 LOCATE ,13: PRINT  WSS : IF PFLG% = 1 THEN LPRINT WSS;" ";
  347. 6280 IF NITER% >= MXITER%                                                               THEN PRINT  MXITER%;" ITERATIONS - PAUSE": GOTO 6920
  348. 6290 IF BRKFLG%                                                                         THEN BRKFLG% = 0 : PRINT "OPERATOR INTERRUPT": GOTO 6920
  349. 6300 COLOR FG2,BG2 : PRINT  "COMPUTING ";: COLOR FG1,BG1
  350. 6305 PRINT " ITERATION #";NITER% + 1;: COLOR FG,BG:PRINT
  351. 6310 '------- DECREASE LAMBDA
  352. 6320 LAMBDA = .31622777# * LAMBDA : INCR = 0
  353. 6340 '------- CALCULATE A=JTJ AND JTF
  354. 6350 FOR I = 1 TO MDI%:A(I) = 0: NEXT
  355. 6360 FOR I = 1 TO NFIT%:GRAD(I) = 0: NEXT
  356. 6370 FOR I = 1 TO NOBS%
  357. 6380     GOSUB 21000
  358. 6390     FMM = (F - DTA(I)) * WGT(I)
  359. 6400     J = 0
  360. 6410     FOR K = 1 TO NPRM%
  361. 6420        IF KFIX(K)=0 THEN J = J+1:DERIV(J) = DERIV(K)*WGT(I)
  362. 6430     NEXT
  363. 6440     FOR J = 1 TO NFIT%
  364. 6450        GRAD(J) = GRAD(J) + DERIV(J) * FMM
  365. 6460        L = J * (J - 1) / 2
  366. 6470        FOR K=1 TO J: A(L+K)=A(L+K)+DERIV(J)*DERIV(K): NEXT
  367. 6480     NEXT
  368. 6490  NEXT
  369. 6500  NDC = NDC + 1
  370. 6510 '------ SAVE "A" MATRIX AND CURRENT PARAMETER VALUES "B"
  371. 6520 FOR I = 1 TO MDI%:C(I) = A(I): NEXT
  372. 6530 FOR I = 1 TO NFIT%:DERIV(I) = B(I): NEXT
  373. 6540 '------ DOCTOR "A" MATRIX DIAGONAL ELEMENTS TO BE:
  374. 6550 '       A = A(1+LAMBDA) + PHI*LAMBDA
  375. 6560 DA = PHI * LAMBDA
  376. 6570 FOR J = 1 TO NFIT%
  377. 6580     G(J) =  - GRAD(J)
  378. 6590     L = J * (J + 1) / 2
  379. 6600     A(L) = C(L) * (1 + LAMBDA) + DA
  380. 6610     K = J - 1
  381. 6620     IF K > 0 THEN  FOR I = 1 TO K:A(L - I) = C(L - I): NEXT
  382. 6630 NEXT
  383. 6650 GOSUB 5000  'Cholesky decomposition
  384. 6660 IF FLG% <> 0 GOTO 6840
  385. 6680 GOSUB 4500   'calculate g (the change in parameters)
  386. 6690 '------ FIND NEW PARAMETERS B = D+G
  387. 6700 '      (N0 COUNTS THE # OF ZERO ELEMENTS IN G)
  388. 6710 N0 = 0 : I = 0
  389. 6720 FOR J = 1 TO NPRM%
  390. 6730     IF KFIX(J) <  > 0 GOTO 6780
  391. 6740     I = I + 1
  392. 6750     B(I) = DERIV(I) + G(I)
  393. 6760     IF  ABS (G(I)) <= ABS (PRCSN*DERIV(I)) THEN N0 = N0 + 1
  394. 6770     PARM(J) = B(I)
  395. 6780  NEXT
  396. 6790  IF N0 = NFIT% GOTO 6890
  397. 6800  OLDWSS = WSS : GOSUB 5500  'update the wss
  398. 6820  IF WSS < OLDWSS GOTO 6250  'next iteration
  399. 6830 '------ last step too big; increase lambda
  400. 6840 IF LAMBDA < PRCSN THEN LAMBDA = PRCSN
  401. 6850 INCR = INCR + 1
  402. 6860 LAMBDA = 10! * LAMBDA
  403. 6870 IF LAMBDA <  = 100000! * LMIN GOTO 6560
  404. 6880 '------ convergence reached
  405. 6890 RI = LAMBDA / LMIN
  406. 6900 CLS : GOSUB 10000 : PRINT  "      CONVERGENCE  REACHED "
  407. 6905 PRINT  "# OF PARAMS FIT = ";NFIT%
  408. 6907 PRINT  "# OF PARAMS CONVERGED = ";N0
  409. 6910 PRINT  "LAMBDA/L(MIN) = ";RI
  410. 6920 NF = NOBS% - NFIT%
  411. 6930 IF NF <  = 0 GOTO 6990
  412. 6940 GOSUB 5700   'calc sum of squares
  413. 6950 SIG =  SQR (SS / NF)
  414. 6960 GOSUB 5900   'calc sum of dta^2
  415. 6970 R =  SQR (SS / R) * 100!
  416. 6980 GOSUB 4000   'compute chi-squared
  417. 6990 IF NITER% <= 0 GOTO 7130
  418. 7000 FOR J = 1 TO MDI%:A(J) = C(J): NEXT : GOSUB 5000
  419. 7010 IF FLG% = 0 GOTO 7060
  420. 7020 FOR J = 1 TO NPRM%
  421. 7030     IF KFIX(J) = 0 THEN PARM(J) = 999.0001
  422. 7040 NEXT
  423. 7050 GOTO 7130
  424. 7060 K = 0
  425. 7070 FOR I4 = 1 TO NPRM%
  426. 7080     IF KFIX(I4) <> 0 GOTO 7120
  427. 7090     K = K + 1
  428. 7100     FOR J = 1 TO NFIT%:G(J) = 0!: NEXT :G(K) = 1!
  429. 7110     GOSUB 4500:DPARM(I4) =  SQR ( ABS (G(K))) * SIG
  430. 7120 NEXT
  431. 7130 RETURN
  432. 9000 'list the current values of the parameters ----------------
  433. 9010 IF NPRM% < 1 THEN  PRINT  "NO PARAM ENTERED": RETURN
  434. 9020 PRINT
  435. 9025 PRINT NPRM%; " FITTING PARAMETERS & THEIR ERRORS" : PRINT
  436. 9030 IF PFLG%= 1 THEN LPRINT : LPRINT NPRM%;                                                       " FITTING PARAMETERS AND THEIR ERRORS":LPRINT
  437. 9040 FOR I = 1 TO NPRM%
  438. 9050     IF PFLG% = 0 GOTO 9080
  439. 9060     LPRINT  I;") ";PARM(I);
  440. 9070     IF KFIX(I) = 0 THEN  LPRINT  " +- ";DPARM(I)                                                   ELSE LPRINT "  (FIXED)  "
  441. 9080     PRINT  I;") ";PARM(I);
  442. 9090     IF KFIX(I) = 0 THEN  PRINT  " +- "; DPARM(I)                                                   ELSE PRINT "  (FIXED)  "
  443. 9100 NEXT
  444. 9110 PRINT : RETURN
  445. 9999 ' ----------- FUNCTION KEY and HELP SUBROUTINES -----------
  446. 10000 'display "help" / title line -----------------------------
  447. 10005 XCURSOR = CSRLIN : YCURSOR = POS(0) : LOCATE 25,1
  448. 10010 IF SHOWTITLE THEN COLOR FG2,BG:PRINT  TITL$;                                    ELSE COLOR FG2,BG : PRINT " 1- SHOW FUNCTION KEYS";                             "   2- SHOW TITLE   3- BREAK   4- QUICK-PLOT";:COLOR FG,BG
  449. 10020 LOCATE XCURSOR,YCURSOR : RETURN
  450. 10200 SHOWTITLE=0: GOSUB 10000: RETURN 'display "HELP" line---F1
  451. 10300 SHOWTITLE = -1: GOSUB 10000: RETURN 'display title -----F2
  452. 10310 BRKFLG% = 1 : RETURN  'user interrupt (break function) -F3
  453. 10400 CLS : PRINT TAB(15); "QUICK PLOT" ' --------------------F4
  454. 10410 PRINT TAB(10); " Y "; YSTYL$; " VS."; " X ";XSTYL$;" PLOT"
  455. 10420 IF PLOTSET% = 0 THEN                                                              PRINT "PLOT NOT DEFINED - FIRST USE MENU OPTION #4" :                          PRINT"PRESS ANY KEY TO CONTINUE":WHILE INKEY$="":WEND:                         RETURN
  456. 10430 GOSUB 16000   'scale data log or linear
  457. 10440 GOSUB 16300   'calculate fitting function
  458. 10450 IF XNEG THEN                                                                       PRINT "WARNING: NEGATIVE X VALUE. X LOG PLOT INVALID":                          XNEG=0:GOTO 10470
  459. 10460 GOSUB 16500 ' PLOT
  460. 10470 SCREEN 0 : RETURN
  461. 15000 'Plot data points and fitting function -------------------
  462. 15010 'Plot data points and fitting function -------------------
  463. 15020 CLS
  464. 15030 IF NPRM% < 1 AND NOBS% < 1 THEN PRINT  TAB(20),                                            "NO DATA POINTS OR PARAMETERS ENTERED":GOTO 400
  465. 15040 FOR I = 1 TO NPRM%
  466. 15050    IF PARM(I) < > 0 THEN GOTO 15080
  467. 15060 NEXT
  468. 15070 PRINT "     WARNING: all parameters are currently = 0"
  469. 15080 RESPLOT = 0   'FLAG to indicate plot for residuals
  470. 15090 IF NOBS% < 1 THEN GOTO 15170
  471. 15100 PRINT TAB(15),"CHOOSE PLOT:":PRINT
  472. 15110 PRINT " 1- PLOT DATA AND FITTING FUNCTION"
  473. 15120 PRINT " 2- PLOT RESIDUALS [ABS(FIT MINUS MEASUREMENT)]"
  474. 15130 PRINT " 3- PLOT RESIDUALS [% FIT MINUS MEASUREMENT]"
  475. 15140 INPUT "ENTER CHOICE:";KMND
  476. 15150 IF KMND = 2 THEN RESPLOT = 1 : PCNT% = 0
  477. 15160 IF KMND = 3 THEN RESPLOT = 1 : PCNT% = 1
  478. 15170 PRINT TAB(15),"CHOOSE GRAPH STYLE:"
  479. 15180 PRINT "1-  Y LINEAR VS. X LINEAR"
  480. 15185 PRINT "2-  Y LINEAR VS. X LOG"
  481. 15190 IF RESPLOT = 1 GOTO 15210
  482. 15200 PRINT "3-  Y LOG VS. X LINEAR":PRINT "4-  Y LOG VS. X LOG"
  483. 15210 OLDSTYL%=GSTYL%:PRINT:INPUT "ENTER CHOICE:";GSTYL%
  484. 15220 IF GSTYL% = 0 THEN 400
  485. 15230 IF GSTYL% <> OLDSTYL% THEN PLOTSET% = 0
  486. 15240 XSTYL$="LINEAR"
  487. 15245 IF (GSTYL%=2 OR GSTYL%=4) THEN XSTYL$="LOG"
  488. 15250 YSTYL$="LINEAR"
  489. 15255 IF (GSTYL%=3 OR GSTYL%=4) THEN YSTYL$="LOG"
  490. 15260 CLS:PRINT TAB(15);" Y ";YSTYL$;" VS.";" X ";XSTYL$;" PLOT"
  491. 15270 IF NOBS% <1 THEN 15460
  492. 15280 ' GO SCALE THE DATA INTO PLOTTING ARRAY
  493. 15290 GOSUB 16000  'scale determine minimum and maximum of DTA
  494. 15300 IF XNEG THEN                                                                       PRINT "WARNING: NEGATIVE X VALUE. X LOG PLOT INVALID":                          XNEG=0:GOTO 15170
  495. 15320 PRINT : PRINT  "DATA POINTS: "
  496. 15330 PRINT  "X VALUES range from: ";XP(1);" to ";XP(NOBS%)
  497. 15340 PRINT  "Y VALUES range from: ";MINDTA;" to ";MAXDTA
  498. 15350 IF PLOTSET% = 0  THEN 15420
  499. 15360 PRINT "Current graph boundaries are:"
  500. 15370 PRINT "       xmin = ";XMIN;"  xmax = ";XMAX
  501. 15380 PRINT "       ymin = ";YMIN;"  ymax = ";YMAX
  502. 15390 PRINT:PRINT                                                                     " Do you wish to change any of the graph boundaries (Y/N)"
  503. 15400 INPUT "(null to leave unchanged)";T1$
  504. 15410 IF T1$ = "" OR T1$ = "N" THEN 15460
  505. 15420 PRINT:INPUT "ENTER GRAPH XMIN,XMAX,YMIN,YMAX ";                                             XMIN,XMAX,YMIN,YMAX
  506. 15425 IF XMIN = 0 AND XMAX = 0 GOTO 420
  507. 15430 PRINT:PRINT "ENTER xinterval,yinterval";                                                           "(THE # OF INTERVALS ON THE x,y AXIS)";
  508. 15435 INPUT XINTERVAL, YINTERVAL
  509. 15440 PLOTSET% = 1
  510. 15450 IF RESPLOT = 1 OR NPRM%=0                                                          THEN XINC = 2*(XMAX-XMIN)/MXOBS%:GOTO 15720
  511. 15460 PRINT:PRINT "PLOT OF FITTED FUNCTION";                                                              "(FROM CURRENT VALUE OF PARAMETERS):"
  512. 15470 IF NOBS% < 1 AND PLOTSET% = 0 GOTO 15530
  513. 15480 PRINT "The function will be plotted in the range ";                                                                      XMIN; " to ";XMAX
  514. 15490 INPUT "ENTER X INCREMENT for plotting the function";XINC
  515. 15500 IF XINC <=0 THEN PRINT "INVALID ENTRY":GOTO 15480
  516. 15510 IF (XMAX-XMIN)/XINC > MXOBS% THEN PRINT                                          "too many points, must be fewer than ";MXOBS%: GOTO 15480
  517. 15520 GOTO 15570
  518. 15530 PRINT "       NO DATA CURRENTLY ENTERED"
  519. 15540 PRINT:PRINT "ENTER XMIN,XMAX AND X-INCREMENT"
  520. 15550 INPUT "for evaluating fitting function";XMIN,XMAX,XINC
  521. 15560 IF (XMAX-XMIN)/XINC > MXOBS% THEN PRINT                                          "too many points, must be fewer than ";MXOBS%: GOTO 15530
  522. 15570 GOSUB 16300 ' evaluate function
  523. 15580 PRINT:PRINT "MINIMUM of fitted function= ";MINYP
  524. 15590 PRINT "MAXIMUM of fitted function= ";MAXYP
  525. 15600 PRINT:PRINT "current graph boundaries are:"
  526. 15610 PRINT "       xmin = ";XMIN;"  xmax = ";XMAX
  527. 15620 PRINT "       ymin = ";YMIN;"  ymax = ";YMAX
  528. 15630 IF NOBS%<1 THEN 15670
  529. 15640 PRINT
  530. 15645 PRINT "do you wish to change THESE GRAPH BOUNDARIES?(Y/N)"
  531. 15650 INPUT "(null to leave unchanged)";T1$
  532. 15660 IF T1$ = "" OR T1$ = "N" THEN 15720
  533. 15670 PRINT:PRINT "ENTER NEW YAXIS BOUNDARY VALUES -- YMIN,YMAX"
  534. 15680 INPUT "(ENTER , TO CHANGE ALL GRAPH BOUNDARIES)";T1$,T2$
  535. 15690 IF T1$ = "" THEN 15420
  536. 15700 YMIN = VAL(T1$):YMAX = VAL(T2$)
  537. 15710 PRINT: PRINT "ENTER xinterval,yinterval";                                                    " (THE # OF INTERVALS ON THE x,y AXIS)"
  538. 15715 INPUT XINTERVAL, YINTERVAL
  539. 15720 GOSUB 16500 'go and plot
  540. 15730 SCREEN 0 : GOTO 420   'finished plotting
  541. 16000 'Scale data & find the minima and maxima of data in YP(I)
  542. 16010 XNEG = 0
  543. 16020 MINDTA = 9.999999E+37: MAXDTA = -9.999999E+37
  544. 16030 FOR I = 1 TO NOBS%
  545. 16040 XP(I) = X(I): YP(I) = DTA(I)
  546. 16050 IF RESPLOT <> 1 THEN 16090
  547. 16060 NORM = DTA(I):IF NORM = 0 THEN NORM = 1E-10
  548. 16070 GOSUB 20000: YP(I) = FUNCTN - DTA(I) 'yes - calc residuals
  549. 16080 IF PCNT% <>0 THEN YP(I) = YP(I)/NORM ' plot % residual?
  550. 16090 IF (GSTYL%=1 OR GSTYL%=3) GOTO 16110 ' xlog plot?
  551. 16100 IF X(I)>0 THEN XP(I) = FNLGT(XP(I)) ELSE XNEG = 1 : RETURN
  552. 16110 IF (GSTYL%=1 OR GSTYL%=2) GOTO 16130 ' ylog plot?
  553. 16120 IF DTA(I)>0 THEN YP(I)=FNLGT(YP(I)) ELSE YP(I) = YMIN
  554. 16130 IF YP(I) >= MAXDTA THEN MAXDTA = YP(I)
  555. 16140 IF YP(I) <= MINDTA THEN MINDTA = YP(I)
  556. 16150 NEXT
  557. 16160 RETURN
  558. 16290 'calculate function using current parm values -----------
  559. 16300 PRINT:PRINT "CALCULATING FUNCTION...."
  560. 16305 IF NPRM% < 1 THEN RETURN
  561. 16310 I=0: KK=0: MINYP = 9.9E+37: MAXYP = -9.9E+37
  562. 16320 FOR XPT=XMIN TO XMAX STEP XINC
  563. 16330     KK=KK+1: FXP(KK) = XPT
  564. 16340     X(I) = XPT :IF (GSTYL%=1 OR GSTYL%=3) GOTO 16360
  565. 16350     X(I)=10^XPT ' FOR X-LOG PLOT
  566. 16360     GOSUB 20000
  567. 16370     FYP(KK) = FUNCTN 'FOR NON-LOG Y PLOT
  568. 16380     IF (GSTYL%=1 OR GSTYL%=2) GOTO 16400
  569. 16390     IF FUNCTN > 0 THEN FYP(KK)=FNLGT(FUNCTN)                                                      ELSE FYP(KK) = YMIN
  570. 16400     IF FYP(KK) >= MAXYP THEN MAXYP = FYP(KK)
  571. 16410     IF FYP(KK) <= MINYP THEN MINYP = FYP(KK)
  572. 16420 NEXT XPT
  573. 16430 KMAX = KK
  574. 16440 RETURN
  575. 16500 'plotting subroutine -- first axes and tick marks --------
  576. 16510 DELTAX = XMAX-XMIN: DELTAY = YMAX-YMIN: LETWID=8: LETHGT=8
  577. 16515 'hi-res screen & choose color 2 -for pc & 100% compatibles
  578. 16520 SCREEN 2:DEF SEG = &H0: OUT &H3D9,2
  579. 16530 LINE (XMIN.PIX,YMIN.PIX) - (XMAX.PIX,YMAX.PIX),,B
  580. 16540 IF XINTERVAL < 1 THEN XINTERVAL = 5
  581. 16550 ' X TICKS
  582. 16560 MINCOL = 0
  583. 16570 FOR TICK = 0 TO XINTERVAL
  584. 16580     TICKPLACE = XMIN.PIX + CINT(DELTAX.PIX*TICK/XINTERVAL)
  585. 16590     LINE(TICKPLACE,YMIN.PIX-2) - (TICKPLACE,YMIN.PIX+2)
  586. 16600     LINE(TICKPLACE,YMAX.PIX-2) - (TICKPLACE,YMAX.PIX+2)
  587. 16610     TICCOL = INT(TICKPLACE/LETWID)-1
  588. 16620     VAL.LAB = XMIN + (TICK/XINTERVAL)*DELTAX
  589. 16630     IF (GSTYL%=2 OR GSTYL%=4) THEN VAL.LAB = 10^VAL.LAB
  590. 16640     TICLAB = VAL.LAB:GOSUB 18000
  591. 16650     TICCOL = TICCOL - LEN(FORMAT$)/2 + 2
  592. 16660     IF TICCOL + LEN(FORMAT$) > 80 THEN 16700
  593. 16670     IF TICCOL < MINCOL GOTO 16700
  594. 16680     LOCATE 25,TICCOL: MINCOL = TICCOL + LEN(FORMAT$)
  595. 16690     PRINT USING FORMAT$;TICLAB;
  596. 16700 NEXT TICK
  597. 16710 IF YINTERVAL < 1 THEN YINTERVAL = 5
  598. 16720 ' Y TICKS
  599. 16730 OLDROW = 26
  600. 16740 FOR TICK = 0 TO YINTERVAL
  601. 16750     TICKPLACE = YMIN.PIX - CINT(DELTAY.PIX*TICK/YINTERVAL)
  602. 16760     LINE(XMIN.PIX-2,TICKPLACE) - (XMIN.PIX+2,TICKPLACE)
  603. 16770     LINE(XMAX.PIX-2,TICKPLACE) - (XMAX.PIX+2,TICKPLACE)
  604. 16780     TICROW = INT (TICKPLACE/LETHGT) + 1
  605. 16790     IF TICROW >= OLDROW-1 THEN GOTO 16860
  606. 16800     VAL.LAB = YMIN + (TICK/YINTERVAL) * DELTAY
  607. 16810     IF (GSTYL%=3 OR GSTYL%=4) THEN VAL.LAB = 10^VAL.LAB
  608. 16820     TICLAB = VAL.LAB :GOSUB 18000
  609. 16830     TICCOL = 8 - LEN(FORMAT$):IF TICCOL < 1 THEN TICCOL=1
  610. 16840     LOCATE TICROW,TICCOL : OLDROW = TICROW
  611. 16850     PRINT USING FORMAT$;TICLAB;
  612. 16860 NEXT TICK
  613. 17000 'plot points and function --------------------------------
  614. 17010 FOR I = 1 TO NOBS%
  615. 17020    CIRCLE (FNSX(XP(I)),FNSY(YP(I))),4
  616. 17030 NEXT
  617. 17040 IF RESPLOT = 1 OR NPRM% < 1 GOTO 17120 ' no function plot
  618. 17050 IWENT%=0
  619. 17060 FOR K = 1 TO KMAX
  620. 17070    XPT = FNSX(FXP(K)): YPT = FNSY(FYP(K))
  621. 17080    IF XPT<0 OR XPT>XTOTAL.PIX OR YPT<0 OR YPT>YTOTAL.PIX                              THEN 17110 ' pt out of range
  622. 17090    IF IWENT% = 0 THEN PSET (XPT,YPT) : IWENT%= 1
  623. 17100    LINE -(XPT,YPT)
  624. 17110 NEXT K
  625. 17120 GTITL$ = TITL$: PCNT$ = "": IF PCNT% = 1 THEN PCNT$="(%)"
  626. 17130 IF RESPLOT = 1 THEN GTITL$ = GTITL$+" - RESIDUALS"+PCNT$
  627. 17140 LOCATE 1,50-LEN(GTITL$):PRINT GTITL$;
  628. 17150 WHILE INKEY$ = "" : WEND : RETURN
  629. 18000 'subroutine to format a label for axis -------------------
  630. 18020 FORMAT$="#"
  631. 18030 IF TICLAB=0 THEN RETURN
  632. 18040 ORDER=INT(LOG(ABS(TICLAB))/2.30258)
  633. 18050 IF ABS(ORDER)>=4 THEN FORMAT$ = "##.#^^^^" : RETURN
  634. 18060 IF ORDER<0 THEN FORMAT$=FORMAT$+".": GOTO 18080
  635. 18070 FORMAT$=FORMAT$+STRING$(ORDER+1,"#")+"."
  636. 18080 T=TICLAB
  637. 18090 TOL=ABS(TICLAB/100000!)
  638. 18100  FP=ABS(T-FIX(T)): IF 1-FP<FP THEN FP=1-FP
  639. 18110 WHILE FP>TOL
  640. 18120   T=T*10: TOL=TOL*10
  641. 18130   FP=ABS(T-FIX(T)): IF 1-FP<FP THEN FP=1-FP
  642. 18140   FORMAT$=FORMAT$+"#"
  643. 18150 WEND
  644. 18160 RETURN
  645. 19000 'error trap ----------------------------------------------
  646. 19010 IF RETCOD = 0 THEN RETCOD = -1 :ERRCNT = 1: GOTO 19030
  647. 19020 ERRCNT = ERRCNT + 1 ' error flagged again in same routine
  648. 19030 IF ERR=5 THEN PRINT "ILLEGAL FUNCTION CALL"
  649. 19040 IF ERR=6 THEN PRINT "OVERFLOW ERROR"
  650. 19050 IF ERR=53 OR ERR=54 THEN PRINT "FILE NOT FOUND"
  651. 19060 IF ERR=61 THEN PRINT "DISK FULL ERROR"
  652. 19070 IF ERR=70 OR ERR=71                                                                THEN PRINT "DISK NOT READY OR WRITE PROTECTED"
  653. 19080 PRINT "ERROR #";ERR;" ON LINE NUMBER ";ERL
  654. 19090 IF ERRCNT > 5 THEN PRINT:PRINT                                                         "ROUTINE TERMINATING DUE TO ERROR COUNT ": GOTO 400
  655. 19100 PRINT "PRESS ANY KEY TO CONTINUE...":WHILE INKEY$="":WEND
  656. 19110 RESUME NEXT
  657. 20000 'subroutine to evaluate the function chosen to be fit
  658. 20005 'The function code here is for the video game problem
  659. 20010 FUNCTN = PARM(1) * (1 - EXP(-PARM(2) * X(I)))
  660. 20020 RETURN
  661. 20050 '
  662. 21000 'subroutine to evaluate the derivative of the function
  663. 21005 'with respect to each parameter. Only 2 parameters here.
  664. 21010 GOSUB 20000
  665. 21020 F = FUNCTN
  666. 21030 DERIV(1) = F / PARM(1)
  667. 21040 DERIV(2) = X(I) * (PARM(1)-F)
  668. 21050 RETURN
  669.